#去除批次效应
setwd("E:\\1.Ovarian_Treg\\1.picture\\3.Batch\\BATCH2")
library(limma)
library(sva)
x=list.files()
x
gse140082=read.table("gse140082.txt",sep = "\t",header = T,check.names = F)
gse140082=gse140082[!duplicated(gse140082$geneNames),]
rownames(gse140082)=gse140082$geneNames
gse140082=gse140082[,-1]
time140082=read.table("time140082.txt",sep = "\t",row.names = 1,header = T,check.names = F)
gse140082_luadexp=gse140082[,row.names(time140082)]
save(gse140082_luadexp,file = "gse140082_luadexp380.rdata")
gse140082_luadexp=cbind(id=row.names(gse140082_luadexp),gse140082_luadexp)
write.table(gse140082_luadexp,file = "gse140082_luadexp380.txt",sep = "\t",row.names = F)


gse14764=read.table("gse14764.txt",sep = "\t",header = T,check.names = F)
gse14764=gse14764[!duplicated(gse14764$geneNames),]
rownames(gse14764)=gse14764$geneNames
gse14764=gse14764[,-1]
time14764=read.table("time14764.txt",sep = "\t",row.names = 1,header = T,check.names = F)
gse14764_luadexp=gse14764[,row.names(time14764)]
save(gse14764_luadexp,file = "gse14764_luadexp80.rdata")
gse14764_luadexp=cbind(id=row.names(gse14764_luadexp),gse14764_luadexp)
write.table(gse14764_luadexp,file = "gse14764_luadexp80.txt",sep = "\t",row.names = F)

gse17260=read.table("gse17260.txt",sep = "\t",header = T,check.names = F)
gse17260=gse17260[!duplicated(gse17260$geneNames),]
rownames(gse17260)=gse17260$geneNames
gse17260=gse17260[,-1]
time17260=read.table("time17260.txt",sep = "\t",row.names = 1,header = T,check.names = F)
gse17260_luadexp=gse17260[,row.names(time17260)]
save(gse17260_luadexp,file = "gse17260_luadexp110.rdata")
gse17260_luadexp=cbind(id=row.names(gse17260_luadexp),gse17260_luadexp)
write.table(gse17260_luadexp,file = "gse17260_luadexp110.txt",sep = "\t",row.names = F)

# gse26193=read.table("gse26193.txt",sep = "\t",header = T,check.names = F)
# gse26193=gse26193[!duplicated(gse26193$geneNames),]
# rownames(gse26193)=gse26193$geneNames
# gse26193=gse26193[,-1]
# time26193=read.table("time26193.txt",sep = "\t",row.names = 1,header = T,check.names = F)
# gse26193_luadexp=gse26193[,row.names(time26193)]
# save(gse26193_luadexp,file = "gse26193_luadexp107.rdata")
# gse26193_luadexp=cbind(id=row.names(gse26193_luadexp),gse26193_luadexp)
# write.table(gse26193_luadexp,file = "gse26193_luadexp107.txt",sep = "\t",row.names = F)

gse32062=read.table("gse32062.txt",sep = "\t",header = T,check.names = F)
gse32062=gse32062[!duplicated(gse32062$id),]
rownames(gse32062)=gse32062$id
gse32062=gse32062[,-1]
time32062=read.table("time32062.txt",sep = "\t",row.names = 1,header = T,check.names = F)
gse32062_luadexp=gse32062[,row.names(time32062)]
save(gse32062_luadexp,file = "gse32062_luadexp260.rdata")
gse32062_luadexp=cbind(id=row.names(gse32062_luadexp),gse32062_luadexp)
write.table(gse32062_luadexp,file = "gse32062_luadexp260.txt",sep = "\t",row.names = F)


# gse49997=read.table("gse49997.txt",sep = "\t",header = T,check.names = F)
# gse49997=gse49997[!duplicated(gse49997$id),]
# rownames(gse49997)=gse49997$id
# gse49997=gse49997[,-1]
# time49997=read.table("time49997.txt",sep = "\t",row.names = 1,header = T,check.names = F)
# gse49997_luadexp=gse49997[,row.names(time49997)]
# save(gse49997_luadexp,file = "gse49997_luadexp194.rdata")
# gse49997_luadexp=cbind(id=row.names(gse49997_luadexp),gse49997_luadexp)
# write.table(gse49997_luadexp,file = "gse49997_luadexp194.txt",sep = "\t",row.names = F)

# gse63885=read.table("gse63885.txt",sep = "\t",header = T,check.names = F)
# gse63885=gse63885[!duplicated(gse63885$geneNames),]
# rownames(gse63885)=gse63885$geneNames
# gse63885=gse63885[,-1]
# time63885=read.table("time63885.txt",sep = "\t",row.names = 1,header = T,check.names = F)
# gse63885_luadexp=gse63885[,row.names(time63885)]
# save(gse63885_luadexp,file = "gse63885_luadexp101.rdata")
# gse63885_luadexp=cbind(id=row.names(gse63885_luadexp),gse63885_luadexp)
# write.table(gse63885_luadexp,file = "gse63885_luadexp101.txt",sep = "\t",row.names = F)


TCGA=read.table("OV_TPM.txt",sep = "\t",row.names = 1,header = T,check.names = F)
group=sapply(strsplit(colnames(TCGA),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
TCGA=TCGA[,group==0]
colnames(TCGA)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", colnames(TCGA))
TCGA=TCGA[,!duplicated(colnames(TCGA))]#
save(TCGA,file = "TCGA_OV422.rdata")
TCGA=cbind(id=row.names(TCGA),TCGA)
write.table(TCGA,file = "TCGA_OV422.txt",sep = "\t",row.names = F)




sam=intersect(rownames(gse140082_luadexp),rownames(gse14764_luadexp))
sam1=intersect(sam,row.names(gse17260_luadexp))
#sam2=intersect(sam1,row.names(gse26193_luadexp))
sam2=intersect(sam1,row.names(gse32062_luadexp))
# sam4=intersect(sam3,row.names(gse49997_luadexp))
#sam4=intersect(sam3,row.names(gse63885_luadexp))
sam5=intersect(sam2,row.names(TCGA))
save(sam5,file = "samegene.rdata")

#果GEO数据没有取log2,会自动对数据取log2
#对基因取交集,分别得到交集基因在TCGA矩阵和GEO矩阵的表达量
gse140082_luadexp=gse140082_luadexp[sam5,]
gse14764_luadexp=gse14764_luadexp[sam5,]
gse17260_luadexp=gse17260_luadexp[sam5,]
#gse26193_luadexp=gse26193_luadexp[sam5,]
gse32062_luadexp=gse32062_luadexp[sam5,]
#gse63885_luadexp=gse63885_luadexp[sam5,]
TCGA=TCGA[sam5,]
##merge
all=cbind(gse140082_luadexp,gse14764_luadexp)
all=cbind(all,gse17260_luadexp)
#all=cbind(all,gse26193_luadexp)
all=cbind(all,gse32062_luadexp)
#all=cbind(all,gse63885_luadexp)
all=cbind(all,TCGA)
#批次矫正

batchType=c(rep(1,ncol(gse140082_luadexp)),
            rep(2,ncol(gse14764_luadexp)),rep(3,ncol(gse17260_luadexp)),
            rep(4,ncol(gse32062_luadexp)),rep(5,ncol(TCGA)))

outTab=ComBat(all, batchType, par.prior=TRUE)

gse140082bat=outTab[,colnames(gse140082_luadexp)]
gse140082bat[gse140082bat<0]=0
save(gse140082bat,file = "gse140082bat.rdata")
gse140082Tab=rbind(ID=colnames(gse140082bat), gse140082bat)
write.table(gse140082Tab, file="gse140082bat.txt", sep="\t", quote=F, col.names=F)

gse14764bat=outTab[,colnames(gse14764_luadexp)]
gse14764bat[gse14764bat<0]=0
save(gse14764bat,file = "gse14764bat.rdata")
gse14764Tab=rbind(ID=colnames(gse14764bat), gse14764bat)
write.table(gse14764Tab, file="gse14764bat.txt", sep="\t", quote=F, col.names=F)

gse17260bat=outTab[,colnames(gse17260_luadexp)]
gse17260bat[gse17260bat<0]=0
save(gse17260bat,file = "gse17260bat.rdata")
gse17260Tab=rbind(ID=colnames(gse17260bat), gse17260bat)
write.table(gse17260Tab, file="gse17260bat.txt", sep="\t", quote=F, col.names=F)

# gse26193bat=outTab[,colnames(gse26193_luadexp)]
# gse26193bat[gse26193bat<0]=0
# save(gse26193bat,file = "gse26193bat.rdata")
# gse26193Tab=rbind(ID=colnames(gse26193bat), gse26193bat)
# write.table(gse26193Tab, file="gse26193bat.txt", sep="\t", quote=F, col.names=F)

gse32062bat=outTab[,colnames(gse32062_luadexp)]
gse32062bat[gse32062bat<0]=0
save(gse32062bat,file = "gse32062bat.rdata")
gse32062Tab=rbind(ID=colnames(gse32062bat), gse32062bat)
write.table(gse32062Tab, file="gse32062bat.txt", sep="\t", quote=F, col.names=F)


# gse63885bat=outTab[,colnames(gse63885_luadexp)]
# gse63885bat[gse63885bat<0]=0
# save(gse63885bat,file = "gse63885bat.rdata")
# gse63885Tab=rbind(ID=colnames(gse63885bat), gse63885bat)
# write.table(gse63885Tab, file="gse63885bat.txt", sep="\t", quote=F, col.names=F)


TCGAbat=outTab[,colnames(TCGA)]
TCGAbat[TCGAbat<0]=0
save(TCGAbat,file = "TCGAbat.rdata")
TCGATab=rbind(ID=colnames(TCGAbat), TCGAbat)
write.table(TCGATab, file="TCGAbat.txt", sep="\t", quote=F, col.names=F)


##对比pca

library("FactoMineR")
library("ggplot2")
library("factoextra")
pca.plot = function(dat,col){
  
  df.pca <- PCA(t(dat), graph = FALSE)
  fviz_pca_ind(df.pca,
               geom.ind = "point", # show points only (nbut not "text")
               pointsize =3,
               pointshape = 21,
               fill.ind = ind, # color by groups
               palette = ggsci::pal_npg("nrc")(10),# c("#00AFBB", "#E7B800", "#FC4E07")
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups",
               title="")+
    theme_bw() +
    theme(text=element_text(size=9,face="plain",color="black"),
          panel.grid = element_blank(),
          axis.title=element_text(size=16,face="plain",color="black"),
          axis.text = element_text(size=14,face="plain",color="black"),
          legend.title = element_text(size=5,face="plain",color="black"),
          legend.text = element_text(size=5,face="plain",color="black"),
          legend.background = element_rect(),
          legend.position="right")
}
# ind即为后续的Legend

x=c(ncol(gse140082bat),ncol(gse14764bat),ncol(gse17260bat),ncol(gse32062bat),ncol(TCGAbat))
ind <- rep(c("GSE140082","GSE14764","GSE17260",
             "GSE32062","TCGA"),x)
pca.plot(all,factor(ind))  ## 处理前
pca.plot(outTab,factor(ind))  ## 处理后
ggsave("Before_batch.pdf",height = 4,width = 5)
ggsave("After_batch.pdf",height = 4,width = 5)
